Índice GeoSES do município de São Paulo

1. Introdução

De forma simplificada, podemos classificar os dados em qualitativos ou quantitativos. Sendo um dado quantitativo, devemos classificá-lo como absoluto ou relativo. As representações por figuras geométricas proporcionais e por densidade de pontos, são usadas para a representação de valores absolutos ou contagens. Quando os dados se referem a índices, taxas, porcentagens, densidades, ou seja, qualquer resultado de uma divisão, são classificados como relativos. Em Cartografia Temática, o conceito envolvido neste tipo de dado não é o de proporcionalidade (como nos dados absolutos), mas o de ordem. A transcrição gráfica da ordem se faz pelas variáveis visuais valor, cor ordenada ou granulação ordenada. Esta representação é a coroplética, a mais frequentemente utilizada em Cartografia Temática.

Vamos utilizar o índice GeoSES no município de São Paulo para inferirmos sobre o contexto socioeconômico da população, por área de ponderação do Censo Demográfico de 2010.

#para limpar a memória
rm(list = ls())

#para ajustar a língua para português usando a função
Sys.setlocale(category = "LC_ALL", 
              locale = "pt_BR.UTF-8")

#para definir o diretório
setwd("C:/Users/Ligia/OneDrive/POS_GRADUACAO_USP/FLG5153_PRINC_CART_TEM2025/AULA7")

#para verificar o diretório
getwd()

#para carregar o pacote 'pacman'
library(pacman)

#para carregar os pacotes necessários
pacman::p_load(
  tidyverse,
  sf,
  colorspace,
  ggplot2,
  ggspatial,
  gridExtra,
  classInt,
  sp)

Vamos abrir o arquivo shapefile das áreas de ponderação do município de São Paulo.

# Para ler o arquivo shapefile
geoses <- st_read("MSP_AREAP_GEOSES.shp")
## Reading layer `MSP_AREAP_GEOSES' from data source 
##   `C:\Users\Ligia\OneDrive\POS_GRADUACAO_USP\FLG5153-PRINCIPIOS DE CARTOGRAFIA TEMATICA 2023\AULA7\MSP_AREAP_GEOSES.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 310 features and 27 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 313394.9 ymin: 7343708 xmax: 360619.8 ymax: 7416189
## Projected CRS: SIRGAS 2000 / UTM zone 23S

Vamos visualizar as áreas de ponderação usando ggplot do pacote ggplot2, definindo a geometria a partir de geom_sf com o contorno na cor preta e o preenchimento em laranja.

#para plotar o mapa das áreas de ponderação
ggplot()+
  geom_sf(data = geoses, colour = "black", fill = "orange")

Mapas coropléticos

Existem decisões importantes quando vamos fazer um mapa coroplético:

  • o número de classes;

  • como fazer a discretização (dividir os intervalos das classes);

  • qual paleta de cores utilizar.

Todos esses aspectos influirão na percepção visual do padrão espacial resultante. Vamos olhar o resumo das variáveis em nosso banco de dados.

O número de classes

Este aspecto afeta diretamente a capacidade de identificar as classes em um mapa. Um mapa pode ser do tipo “sem classes” ou do tipo “com classes”. Os mapas sem classes atribuem uma tonalidade a cada valor do conjunto de dados. É muito interessante para termos ideia do padrão geral do mapa. No entanto, o leitor dificilmente conseguirá identificar o valor de cada tonalidade porque a legenda será apresentada em coluna única. Assim, os mapas divididos em classes são mais utilizados.

Ao definir o número de classes, se utilizarmos poucas classes (menos de 4), podemos perder as situações intermediárias do fenômeno representado, resultando em grande simplificação. Se usarmos mais de 9 classes, por outro lado, o leitor terá dificuldade de identificar as tonalidades em uma escala monocromática. O intervalo entre 4 e 7 classes tem sido apontado como o mais indicado. Podemos utilizar a fórmula de Sturges para orientação sobre o número de classes (k) (CAUVIN et al., 2007). Na fórmula a seguir, n significa o número de unidades geográficas (no nosso caso, o número de áreas de ponderação = 310).

\[ k = 1 + 3,3 ∗ log(n) \] Como vimos com a função ´summary´, o shapefile de ‘geoses’ tem 310 áreas de ponderação. Aplicando a fórmula de Sturges, teríamos 9,22 classes. De acordo com a fórmula, com 70 unidades geográficas já atingiríamos o limite de 7 classes. Assim, como regra geral, acima de 70 unidades, podemos usar 7 classes, para uma paleta monocromática. Se a paleta tiver duas cores opostas (e.g., azul e vermelho), poderíamos ter até 7 classes em tons de azul e até sete em tons de vermelho.

Divisão dos intervalos de classes

A segunda decisão necessária diz respeito a como dividir os intervalos das classes. Para isso, precisaremos recorrer a conceitos muito básicos de estatística descritiva. Na prática, precisamos avaliar a distribuição dos dados da variável de interesse. A melhor forma de fazer isso é examinarmos o histograma da variável. Um histograma mostra a frequência dos intervalos dos dados e permite avaliar se a distribuição se assemelha a uma distribuição normal (Gaussiana). A distribuição normal é uma curva simétrica em torno do seu ponto médio, apresentando um formato de sino. Quando os dados apresentam simetria, a forma de discretização afetará pouco o padrão espacial do mapa. No entanto, quando há grande assimetria nos dados é necessário um cuidado maior para que a representação cartográfica comunique o padrão espacial do fenômeno de forma mais adequada.

#para gerar um histograma da variável 'GeoSES'
histo_geoses <- ggplot(geoses, aes(x=GeoSES))+ 
  geom_histogram(bins=12, aes(y = after_stat(density)), 
                              colour = "blue", 
                              fill="lightblue", 
                              na.rm = T)+
  geom_vline(xintercept = mean(geoses$GeoSES, na.rm = T), 
             linewidth = 1, 
             colour = "orange", 
             linetype = "dashed")+
  stat_function(fun = dnorm, 
                colour = "red", 
                linewidth = 1,
                args = list(mean = mean(geoses$GeoSES, na.rm = T),
                            sd = sd(geoses$GeoSES, na.rm = T))
  )
#para ver o histograma
histo_geoses

Vemos que existem mais áreas de ponderação com valores de GeoSES abaixo de zero.

Vamos elaborar o mapa sem classes.

#para elaborar o mapa sem classes
mapa_semclasses <- ggplot()+
  geom_sf(data = geoses, 
          mapping = aes(fill = GeoSES)) +
  scale_fill_continuous_diverging("Red-Green") +
  ggspatial::annotation_scale(location = "br", width_hint = 0.2, 
                              line_width = 0.5, height = unit(0.1,"cm"))

#para visualizar o mapa
mapa_semclasses

Vamos exportar o mapa elaborado para posterior comparação com o mapa com classes.

#para exportar o mapa como figura jpg
ggsave("mapa_semclasses.jpg", 
       width = 8, 
       height = 8, 
       dpi = 300)

O pacote classInt apresenta a função jenks.test que avalia a acurácia dos intervalos criados segundo os diferentes algoritmos.

Vamos calcular objetos do tipo classInt para avaliar a acurácia dos intervalos para cada tipo de discretização. Existem 2 índices: índice de Jenks e o Tabular Accuracy Index (TAI).

Vamos criar objetos do tipo classIntervals para GeoSES para cada estilo disponível.

#para criar os intervalos para as discretizações disponíveis
geoses_fixed <- classIntervals(geoses$GeoSES, n=7, style="fixed",
 fixedBreaks=c(-1, -0.66, -0.33, 0.0, 0.33, 0.66, 1))

geoses_sd <- classIntervals(geoses$GeoSES, n=7, style="sd")

geoses_equal <- classIntervals(geoses$GeoSES, n=7, style="equal")

geoses_quant <- classIntervals(geoses$GeoSES, n=7, style="quantile")

geoses_kmeans <- classIntervals(geoses$GeoSES, n=7, style="kmeans")

geoses_pretty <- classIntervals(geoses$GeoSES, n=7, style="pretty")

geoses_fisher <- classIntervals(geoses$GeoSES, n=7, style="fisher")

geoses_jenks <- classIntervals(geoses$GeoSES, n=7, style="jenks")
#para rodar e visualizar o teste para cada discretização
print(jenks.tests(geoses_fixed))
##        # classes  Goodness of fit Tabular accuracy 
##         6.000000         0.953987         0.773591
print(jenks.tests(geoses_sd))
##        # classes  Goodness of fit Tabular accuracy 
##       10.0000000        0.9802571        0.8518538
print(jenks.tests(geoses_equal))
##        # classes  Goodness of fit Tabular accuracy 
##        7.0000000        0.9633570        0.7984705
print(jenks.tests(geoses_kmeans))
##        # classes  Goodness of fit Tabular accuracy 
##        7.0000000        0.9668289        0.8201267
print(jenks.tests(geoses_pretty))
##        # classes  Goodness of fit Tabular accuracy 
##        4.0000000        0.9020182        0.6771809
print(jenks.tests(geoses_fisher))
##        # classes  Goodness of fit Tabular accuracy 
##        7.0000000        0.9743273        0.8384473
print(jenks.tests(geoses_jenks))
##        # classes  Goodness of fit Tabular accuracy 
##        7.0000000        0.9743273        0.8384473

Os resultados mostram que o desvio-padrão (geoses_sd) com 10 classes apresenta a melhor acurácia (GOF de 0.98 e TAI de 0.85). Vamos visualizar os intervalos.

#para ver os intervalos de 'geoses_sd'
geoses_sd$brks
##  [1] -1.01434275 -0.79970340 -0.58506404 -0.37042469 -0.15578533  0.05885402
##  [7]  0.27349338  0.48813273  0.70277209  0.91741145  1.13205080

Vamos criar uma coluna com o nome ‘geoses_class’.

#para criar uma coluna chamada 'geoses_class' com os valores das quebras de desvio-padrão
geoses <- geoses %>%
  mutate(geoses_class = cut(GeoSES, geoses_sd$brks, include.lowest = TRUE))

Vamos criar um vetor com os valores de intervalos definidos pela função ‘classInt’, usando os valores reais do início e fim dos valores (-1 e 1).

#vamos usar os valores de intervalos definidos pela função 'classInt', usando os valores reais do início e fim dos valores (-1 e 1)
legend_values <- c("-1.00 a -0.79", "-0.79 a -0.58", "-0.58 a -0.37", "-0.37 a -0.15", "-0.15 a 0.05",  "0.05 a 0.27", "0.27 a 0.48", "0.48 a 0.70", "0.70 a 0.91", "0.91 a 1.00")

Vamos elaborar o mapa com classes.

#para elaborar o mapa
mapa <- ggplot()+
  geom_sf(data = geoses, 
          aes(fill = geoses_class)) +
  scale_fill_discrete_diverging("Red-Green", 
                                labels = legend_values) +
  labs(fill = "GeoSES") +
  ggspatial::annotation_scale(location = "br", width_hint = 0.2, 
                              line_width = 0.5, height = unit(0.1,"cm"))

#para visualizar o mapa
mapa

#para exportar o mapa como figura jpg
ggsave("mapa.jpg", 
       width = 8, 
       height = 8, 
       dpi = 300)

Vamos comparar o resultado visual dos dois mapas.

#para colocar os dois mapas lado a lado
mapa2 <- grid.arrange(mapa_semclasses, 
             mapa, 
             ncol=2
)

Vamos exportar os 2 mapas juntos.

#para exportar o mapa como figura jpg
ggsave("mapa2.jpg", 
       mapa2,
       device = "jpeg",
       width = 8, 
       height = 3, 
       dpi = 300)

Análise do padrão espacial do GeoSES

Vamos calcular o I de Moran para o GeoSES para avaliar a tendência de agrupamento global no município de São Paulo.

#para carregar os pacotes
library(devtools)
library(spdep)

Cálculo do coeficiente I de Moran

Esta metodologia para cálculo do coeficiente I de Moran segue a proposta por https://mgimond.github.io/es214_support_tutorials/moranI/Mapping_and_Morans.html.

Passo 1: Definir os polígonos vizinhos

Vamos definir uma matriz de vizinhança do tipo rainha (“queen”).

#para criar a matriz de vizinhança do tipo 'queen'
nb <- poly2nb(geoses, queen=TRUE)

Passo 2: Atribuir pesos aos vizinhos

#para atribuir os pesos
lw <- nb2listw(nb, style="W", zero.policy=TRUE)

#para ver os pesos dos vizinhos do polígono 5
lw$weights[5]
## [[1]]
## [1] 0.25 0.25 0.25 0.25

Passo 3: Computar os valores do GeoSES

#para calcular a defasagem espacial
geoses.lag <- lag.listw(lw, geoses$GeoSES)
#para plotar os índices e os índices ponderados pelos vizinhos
plot(geoses.lag ~ geoses$GeoSES, pch=16, asp=1)
abline(lm(geoses.lag ~ geoses$GeoSES), col="blue")

Passo 4: Computar a estatística I de Moran (global)

#para calcular o coeficiente I de Moran
I <- moran(geoses$GeoSES, lw, length(nb), Szero(lw))[1]
#para ver o valor
I
## $I
## [1] 0.845167

Passo 5: Realizando um teste de hipóteses

A hipótese a ser testada é a de que “os valores dos índices são distribuídos aleatoriamente ao longo do território seguindo um processo completamente aleatório”. Há dois métodos para testar esta hipótese: um analítico e um por meio de simulação do tipo Monte Carlo.

Método analítico

Para rodar a análise usando o método analítico, usaremos a função moran.test.

#para rodar o teste
moran.test(geoses$GeoSES,lw, alternative="greater")
## 
##  Moran I test under randomisation
## 
## data:  geoses$GeoSES  
## weights: lw    
## 
## Moran I statistic standard deviate = 24.305, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.845167046      -0.003236246       0.001218484

Método de simulação Monte Carlo

A vantagem do método analítico é que ele é rápido. No entanto, pode ser sensível a polígonos de distribuição irregular. Uma abordagem mais segura é rodar uma simulação do tipo Monte Carlo usando a função moran.mc. O número de simulações é definido pelo parâmetro ‘nsim’. Vamos usar valores de 999 vezes.

A função moran.mc tem outro parâmetro chamado ‘alternative =’. Este parâmetro tem três valores possíveis: “greater” (the default), “less” e “two.sided”. A escolha será ditada pelo lado da distribuição que nós queremos computar para o valor de p. Se o valor observado de Moran estiver do lado direito da distribuição esperada, nós vamos querer adotar a opção “greater” que focará na cauda superior da distribuição. Se estiver do lado esquerdo, nós vamos querer escolher “less” para focar na cauda inferior da distribuição. No nosso caso, o valor I de Moran é 0.85, positivo, então vamos escolher “greater” para o parâmetro.

#para rodar a simulação de Monte Carlo para o cálculo de I de Moran
MC<- moran.mc(geoses$GeoSES, lw, nsim = 999, alternative = "greater")

# para ver os resultados
MC
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  geoses$GeoSES 
## weights: lw  
## number of simulations + 1: 1000 
## 
## statistic = 0.84517, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater
# plotar a distribuição nula (note que este é um gráfico de densidade e não um histograma)
plot(MC)

A curva mostra a distribuição dos valores que poderíamos esperar se a taxa fosse distribuída aleatoriamente pelas áreas de ponderação. Note que o resultado de 0.85 cai fora no lado direito da distribuição sugerindo que os valores dos índices estão agrupadas (‘clustered’). Um valor positivo de I de Moran sugere ‘clustering’ enquanto um valor negativo sugere dispersão.

# para plotar o diagrama de autocorrelação de Moran
moran.plot(geoses$GeoSES, lw)

A simulação de Monte Carlo gerou um valor de p de 0.001, sugerindo que haveria 0.001 % de chance de estarmos errados em rejeitar a hipótese nula ou que há 0,1% de chance que o padrão observado é consistente com um processo aleatório. O mesmo conjunto de dados pode apresentar resultados de p diferentes, já que resulta de simulações. Quanto maior o número de simulações, mais estável fica o resultado.

Podemos concluir, portanto, que o índice GeoSES apresenta forte padrão espacial global significativo.

Cálculo do Moran Local (LISA - “Local Indicator of Spatial Association”)

Vamos usar o pacote ‘rgeoda’ para a análise de agrupamentos locais. Veja este tutorial em https://geodacenter.github.io/rgeoda/articles/rgeoda_tutorial.html.

#para carregar o pacote
library(rgeoda)
## Loading required package: digest
## 
## Attaching package: 'rgeoda'
## The following object is masked from 'package:spdep':
## 
##     skater
#para calcular uma matriz de pesos do tipo 'queen'
queen_w <- queen_weights(geoses)

#para selecionar a coluna GeoSES
geo = geoses["GeoSES"]

#para calcular o moran local do GeoSES
lisa_geo <- local_moran(queen_w, geo)

#para definir as cores
lisa_colors <- lisa_colors(lisa_geo)

#para definir os rótulos 
lisa_labels <- lisa_labels(lisa_geo)

#para definir os clusters
lisa_clusters <- lisa_clusters(lisa_geo)

#para plotar o mapa do Moran local
plot(st_geometry(geo), 
     col=sapply(lisa_clusters, function(x){return(lisa_colors[[x+1]])}), 
     border = "#333333", lwd=0.2)
title(main = "Moran Local do GeoSES")
legend('bottomleft', legend = lisa_labels, fill = lisa_colors, border = "#eeeeee")

O mapa LISA nos mostra que existe uma área com altos valores de GeoSES em vermelho, que formam um agrupamento espacial e outras áreas em azul, com baixos valores de GeoSES que se localizam na periferia. Os agrupamentos são significativos estatisticamente. Ou seja, não acontecem por acaso. As áreas em cinza são não significativas, pois estão dentro da média dos valores da cidade. Podemos constatar que existe uma forte segregação residencial em relação ao contexto socioeconômico.

Vamos visualizar o mesmo resultado usando a base do OpenStreetMap.

#para adicionar os clusters na tabela
geoses$moran_cluster <- lisa_clusters

#para converter para character
geoses$moran_cluster <- as.character(geoses$moran_cluster)

# Paleta de cores personalizada #para definir paletas, veja https://r-charts.com/color-palettes/

cores_lisa <- c("white", "red", "blue", "lightblue", "lightpink")

#para carregar os pacotes
library(ggplot2)

# para criar um fator com níveis personalizados e rótulos
geoses$moran_cluster <- factor(geoses$moran_cluster, 
                               levels = c("0", "1", "2", "3", "4"), 
                               labels = c("Não sign.", "Alto-Alto", "Baixo-Baixo", "Baixo-Alto", "Alto-Baixo"))
#para visualizar o mapa de clusters no OpenStreetMap
mapview::mapview(geoses, 
                 zcol = "moran_cluster", 
                 legend = TRUE, 
                 legend.style = "fill", 
                 layer.name = "Índice GeoSES", 
                 col.regions = cores_lisa)

Referências

Appelhans T, Detsch F, Reudenbach C, Woellauer S (2023). mapview: Interactive Viewing of Spatial Data in R. R package version 2.11.2, https://CRAN.R-project.org/package=mapview.

Auguie B (2017). gridExtra: Miscellaneous Functions for “Grid” Graphics. R package version 2.3, https://CRAN.R-project.org/package=gridExtra.

Barrozo, L. V. et al. GeoSES: A socioeconomic index for health and social research in Brazil. PLOS ONE, v. 15, n. 4, p. e0232074, 2020.

Bivand R (2023). classInt: Choose Univariate Class Intervals. R package version 0.4-10, https://CRAN.R-project.org/package=classInt.

Bivand, R. (2022). R Packages for Analyzing Spatial Data: A Comparative Case Study with Areal Data. Geographical Analysis, 54(3), 488-518. doi:10.1111/gean.12319 https://doi.org/10.1111/gean.12319.

Bivand R, Pebesma E, Gomez-Rubio V (2013). Applied spatial data analysis with R, Second edition. Springer, NY. https://asdar-book.org/.

Cauvin, C. Escobar, F., Serradj, A. Cartographie thématique: une nouvelle démarche (Vol. 1). 2007.

Dunnington D (2023). ggspatial: Spatial Data Framework for ggplot2. R package version 1.1.9, https://CRAN.R-project.org/package=ggspatial.

Li, X., Anselin, L. (2023). rgeoda: R Library for Spatial Data Analysis. R package version 0.0.10-4, https://CRAN.R-project.org/package=rgeoda.

Pebesma E, Bivand R (2005). Classes and methods for spatial data in R. R News, 5(2), 9-13. https://CRAN.R-project.org/doc/Rnews/.

Pebesma, E., & Bivand, R. (2023). Spatial Data Science: With Applications in R. Chapman and Hall/CRC. https://doi.org/10.1201/9780429459016

Pebesma, E., 2018. Simple Features for R: Standardized Support for Spatial Vector Data. The R Journal 10 (1), 439-446, https://doi.org/10.32614/RJ-2018-009

Rinker, T. W. & Kurkiewicz, D. (2017). pacman: Package Management for R. version 0.5.0. Buffalo, New York. http://github.com/trinker/pacman

R Core Team (2023). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.

Wickham, H. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2016.

Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). Welcome to the tidyverse. Journal of Open Source Software, 4(43), 1686. doi:10.21105/joss.01686 https://doi.org/10.21105/joss.01686.

Wickham H, Hester J, Chang W, Bryan J (2022). devtools: Tools to Make Developing R Packages Easier. R package version 2.4.5, https://CRAN.R-project.org/package=devtools.

Zeileis A, Fisher JC, Hornik K, Ihaka R, McWhite CD, Murrell P, Stauffer R, Wilke CO (2020). colorspace: A Toolbox for Manipulating and Assessing Colors and Palettes. Journal of Statistical Software, 96(1), 1-49. doi:10.18637/jss.v096.i01 https://doi.org/10.18637/jss.v096.i01.

Zeileis A, Hornik K, Murrell P (2009). Escaping RGBland: Selecting Colors for Statistical Graphics. Computational Statistics & Data Analysis 53(9), 3259-3270. doi:10.1016/j.csda.2008.11.033 https://doi.org/10.1016/j.csda.2008.11.033.